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(57) Abstract: A system and method of creating a filter for 
use with locally dense seismic data is disclosed. The method 
includes obtaining survey geometry characteristics from a lo- 
cally dense seismic survey. A filter is designed which uses 
spatial derivatives of the wavefield of order between (1) and 
the maximum order of spatial derivatives of the wavefield that 
can be estimated within a group. The filter can be designed so 
as to separate up/down going components, p/s components, 
or both up/down and p/s components. Partial derivatives in 
space and time of the wavefield can be calculated, using, for 
example, a taylor series expansion as an approximation. The 
seismic data is filtered by combining estimated near surface 
material properties, the seismic data, and the calculated par- 
tial derivatives. 
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System and Method for Seismic Wavefield Separation 

FIELD OF THE INVENTION : 

The present invention relates to the field of 
5 seismic data acquisition and processing. In 

particular, the invention relates to a system and 
method for seismic wavefield separation. 

BACKGROUND OF THE INVENTION : 

10 Optimal processing, analysis and 

interpretation of land seismic data ideally require 
full information about the wavefield so that the 
wavefield can be separated into its up- and down-going 
and P- and S-components as well as determining phase 

15 and polarity. For 3C acquisition of land surface 

seismic data it is common practice simply to interpret 
the vertical component as the P-section and the 
horizontal components as SV- and SH-sections. This 
* traditional" P/S interpretation is exact for vertical 

20 arrivals. However, as energy is incident away from 
normal incidence angles, this approximation breaks 
down, both because of projections on to all components, 
but also because of a non-unity reflection coefficients 
and mode-conversions at the free-surface. 

25 Exact analytical filter expressions for 

wavefield separation have previously been derived by 
for instance Dankbaar, J. W. M. , 1985, Separation of P- 
and S-waves: Geophys. Prosp., 33, 970 - 986, and these 
have been applied to seismic data in conventional 

30 recording geometries. Unfortunately, statics problems 
severely limit the practical use of these wave- equation 
based techniques . 



- WO 01/53854 PCT/GB01/00186 

- 2 - 

SUMMARY OF THE I3STVENTION : 

Thus, it is an object of the present 
invention to provide a filtering technique that lends 
itself; to be applied within local densely deployed 
5 single sensor receiver groups where statics are 
substantially constant. 

It is a further object of the present 
invention to provide a filtering technique that can be 
implemented efficiently directly in the spatial domain. 

10 According to the invention a method of 

creating a filter for use with locally dense seismic 
data is provided. The method includes obtaining survey 
geometry characteristics from a locally dense seismic 
survey designed to record characteristics of an elastic 

15 or acoustic wavefield. The seismic survey is made up 
of a number of groups of receivers, with each group 
comprising at least three receivers densely spaced from 
each other. According to the method, a filter is 
designed which uses spatial derivatives of the 

20 wavefield of order between 1 and the maximum order of 
spatial derivatives of the wavefield that can be 
estimated within a group. The filter is designed such 
that when combined with data from within a single 
group, the filter separates components of some or all 

25 of the wavefield arriving at the single group. 

The filter can be designed so as to separate 
up/down going components, p/s components, or both 
up/down and p/s components. 

The densely spaced receivers in the group are 

30 preferably spaced apart such that statics in the 

portion of the wavefield of interest are substantially 
constant. More preferably, each of the densely spaced 
receivers in the group are spaced apart by about 2 
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meters or less, or by a distance of about one fifth the 
shortest wavelength of interest or less. 

Partial derivatives of the wavefield are also 
preferably calculated, and this can be done. using a 
5 taylor series expansion as an approximation. According 
to the invention, the seismic data is preferably 
filtered by combining estimated near surface material 
properties, the seismic data, and the calculated 
partial derivatives (both in space and time) . 

10 The filter can also be used to separate 

surface waves or airwave induced ground motion from the 
seismic data. The free surface condition can be used 
to convert vertical derivatives of the wavefield to 
horizontal derivatives of the wavefield. 

15 According to the invention, the seismic 

survey is performed primarily for the purpose 
hydrocarbon reservoir exploration, evaluation or 
characterisation, although other uses can be made. 

The invention is applicable where the near 

20 surface velocities are isotropic or anisotropic. 

BRIEF DESCRIPTION OF THE DRAWINGS : 

Figures la-c show various illustrations of 
25 the elastodynamic representation theorem, according to 

the invention; 

Figure 2 shows a plot of the normalised real 

part of the filter for removing free surface effects 

from divergence in equation (36) ; 
30 Figures 3a and 3b show an example of a 

seismic survey having locally dense receiver groups, 

according to the invention; 
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Figures 4a and 4b show the divergence as 
calculated directly by numerical methods ; 

Figures 5a-d illustrate sections obtained 
. from V3 measurements only; 
5 Figures 6a-d show the results of using a 

wave field separation technique according to a preferred 
embodiment of the invention; 

Figure 7 is a schematic illustration of a 
seismic data acquisition and processing system, 
10 according a preferred embodiment of the invention; and 
Figure 8 shows steps involved in the 
filtering technique, according to preferred embodiments 
of the invention. 

15 

DETAILED DESCRIPTION OF THE INVENTION : 

According to the invention a new approach is 
provided for P/S separation of land surface seismic 
data and for removing the effects of the Earth's free 

20 surface (i.e. up/down separation) . By converting 

vertical spatial derivatives to horizontal derivatives 
using the free- surface condition, the methodology can 
make use of locally dense measurements of the wavefield 
at the free surface to calculate all spatial 

25 derivatives of the wavefield. These can in turn be 

used to compute divergence (P-waves) and curl (S-waves) 
of the wavefield at the free surface. 

The effects of the free surface are 
preferably removed through an up/down separation step 

30 using the elastodynamic representation theorem. This 
results in expressions where spatial filters are 
convolved with recorded data. The filters can be 
successfully approximated so that they fit the local 
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dense acquisition pattern used for the P/S separation 
step. In particular; the simplest approximate 
expression for up-^going P-waves consists of two terms. 
The first term preferably corresponds to the divergence 
5 in the presence of the free surface scaled by a 

material constant. The second term preferably is a 
time derivative of the recorded vertical component 
scaled by a material constant. Hence, a correction is 
added to the "traditional" P- interpretation through the 

10 first term, which improves accuracy for incidence 
angles outside normal incidence. 

In UK Patent Application entitled "System and 
Method for Estimating Seismic Material Properties" (UK 
Patent Application No. 0003410.8) filed on 15 February 

15 2000, incorporated herein by reference, a method is 

provided than can use the volumetric recordings of the 
wavefieid to invert for P and S velocities in the Earth 
in the neighbourhood of a small, closely-spaced array 
of receivers. "Volumetric recording" refers to an 

20 array that approximately encloses a volume of the 
Earth. The quantities estimated are the effective 
velocities of the P- and S-components of the wavefieid 
at any point in time. Hence, these will vary with both 
wave type and wavelength. If estimated for the near 

25 surface Earth structure, such velocities may be useful 
for statics estimation, or for separation of the 
wavefieid into up- and down-going components as 
described herein. 

Implications of the free surface condition 

30 for land seismic recording will be discussed first with 
reference to the fully anisotropic case; The elastic 
constitutive relation relates components of the stress 
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tensor ay in a source free region. to the strain tensor 
components 



(1) 



where ajja are the elastic stiffnesses. Index 
values 1 and 2 correspond to horizontal coordinates xj 
and *2 whereas index value 3 corresponds to the 
vertical direction downwards X3. Using the Voigt 
10 notation, equation (1) can be written as: 
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where C is the symmetric stiffness matrix 
15 with 21 independent components: 
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(3) 



The strain % is related to particle velocity 



20 Vf through 



eij= ^(djVi+diVj), 



(4) 
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where 9i denotes spatial derivatives in the 
xj, X2 or X3 directions, and the dot denotes a time 
derivative. In land acquisition, spatial distributions 
5 of 3C receivers at the surface allow us to compute 
horizontal spatial derivatives of particle velocities 
(or time-derivatives thereof if particle acceleration 
is recorded, etc.). The only information missing in . 
order to know the wavefield completely is the vertical 
10 derivatives of the recorded wavefield. 

The free-surface condition at the Earth's 
surface gives us three additional constraints: 

^3 = 0, (5) 

15 

which are sufficient to allow us to compute 
the remaining velocities given that we make some 
independent estimate of the relevant elastic 
stiffnesses. 

20 In addition, constraints on the relation 

between individual elements of the stress and strain 
tensors could possibly be used to correct for coupling 
or to compute near-surface properties. 

The case of isotropic material properties in 

25 the near-surface environment is of particular interest 
in the land surface seismic case. Using the Voigt 
notation, the stiffness matrix takes the following 
form: 
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(6) 



where A and fi are the Lame constants . 
The constraints imposed by the free-surface 
5 condition (5) become: 



03t>3 = 



(7) 

(8) 
(9) 



where the horizontal derivatives on- the right 
10 hand side are known from the surface measurements . 
Note that the material properties only occur in 
equation (7) and not in equations (8) and (9) . 

Divergence and curl of a wavefield at a free 
surface overlaying a homogeneous isotropic half space 
15 will now be discussed. . This discussion is in the 
context of an isotropic media. The elastic wave 
equation for particle displacement u can be written as: 



20 



pu = f + (A + 2/i)V(V-u)~ juV x (Vx u), 



(10) 



where p is the density and f denotes a 
distribution of body forces. Lamp's theorem states 
that there exist potentials $ and Y of u with the 
following properties: 
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U = V$ + Vx$, (U) 

V * = 0 , (12) 

$ - i + clV 2 $, (13) 
P 

* = t+dV** 3 (14) 



where, c a and cp are the P- and S-velocities, and <j> and yr 
5 are potential components associated with the body 

force. An elastic wavefield u can thus be decomposed 
into its P~ and S-wave components, V<£ and Vx*F, 
respectively. Equations (11) and (12) yield: 

V-u = V 2 $, (15) 
io Vxu = VxVxf (16) 

By measuring the curl and the divergence of 
an elastic wavefield we can thus measure the P- and S- 
wave components separately. 

15 An acquisition pattern comprising tetrahedra 

of 3C measurments can be used to achieve the separation 
of the wavefield into its curl- and divergence- free 
components. See, e.g. Robertsson, J.O.A., and Muyzert, 
E . / 1999, Wavefield separation using a volume 

20 distribution of three component recordings: Geophys. 

Res. Lett., vol. 26, 2821-2824, incorporated herein by 
reference. Such acquisition patterns are described in 
further detail in UK Patent Application No. 9921816 • 6, 
incorporated herein by reference. All spatial 

25 derivatives of the wavefield components can be 
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calculated. Therefore, divergence and curl can be 
calculated explicitly from surface measurements only. 
Equations (7), (8) and (9) give us the following 
expressions for divergence and curl of particle 



5 velocity at a free surface: 

(Vxv)i = 2d 2 v z , (18) 
( V x v) 2 = - 2d 1 v l , (19) 
(V x v) 3 = div 2 - d 2 v 1 . (20) 

First we note that the ratio — - — = 2 (cp/c a ) 



(may be frequency dependent) scales the expression for 

10 the divergence in equation (17) . 

Second, equations (17), (18), (19) and (20) 
contain both the up-going and down-going parts of the 
wavefield. This includes mode conversions at the free 
surface. For instance, the divergence given by (17) 

15 contains not only the desired up-going P-waves, but 
also, the down-going P-to-P reflection, down-going S- 
to-P conversions, etc. Moreover, a plane P-wave which 
is vertically incident on the free-surface will have 
zero divergence (up- and down-going parts interfere 

20 destructively) . Removing the effects of the free 

surface is therefore regarded as an important step in 
the preferred P/S separation technique. 

A technique for Up/down wavefield separation 
using the elastodynamic representation theorem will now 

25 be discussed. The elastodynamic representation theorem 
or Betti's relation can be derived from the equation of 
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20 



motion and the elastic constitutive relations using 
Gauss' theorem. Suppose that we have a volume V 
enclosed by a surface S, and that we wish to calculate 
the displacement of a wavefield u at a point x' in V. ■". 
The displacement u(x') is directly related to the stress 
a and displacement along S, sources in V, as well as 
the displacement Green's tensor Gy and the stress 
Green's tensor functions between S or source points .... 
and x' : 

«»M = £(« i (x)Gi„(x,x')-".(=')«Av„(x,x'))«l5(x) (21) 



where n is the normal unit vector to S, ti(x) 
is the traction across S at point x and fi is the i m 
component of the force from sources within V. 
15 Equation (21) shows the frequency domain expression of 
the representation theorem. In the time domain we 
would also have convolutions in time. The traction t 
along S is defined as 



t = n • (7. (22) 



Now, suppose that the volume V consists of a 
homogeneous elastic medium. Figures la-c show various 
illustrations of the elastodynamic representation 
25 theorem, according to the invention. First, consider 
the closed surface Si +S 2 illustrated in Figure la. In 
this case we assume that we have a source located 
outside V so that it is recorded along Si as an up- 
going mode. Sommerf eldt 1 s radiation condition implies 
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that the contribution from S2 vanishes at infinity. 
Hence, the closed surface integral in equation (21) can 
be replaced by the open surface integral: 

u n (x') = / (ti(x)G lM (x,x') - Ji) (x)« i S ijn (x J x , ))a5( X ) . (23) 



The expression in equation (23) represents an 
up-going wavefield at x' since the field is up-going at 
10 Si. 

Next consider the situation depicted in 
Figure lb. If we in a similar way put the source above 
Si outside V, we obtain an analogous expression to 
equation (23) for the down-going wavefield. 

15 Finally, we consider the situation shown in 

Figure lc. This time we assume that Si follows the 
Earth's surface topography, that no sources (primary or 
secondary) are located in V, and that x' is located 
inf initesimally below the land surface. Again, the 

20 left and right edges of the surface S are at infinity 
and do not contribute to the integral. The 
representation theorem therefore becomes: 



t>n(x') = / (Ho)t<x)G,- n ( X) x') - t/;(x)n x')) dS(x) . (24) 

25 

In equation (24) we have used particle 
velocities V/ instead of displacements Ui since this will 
be more convenient to work with. From the above 
30 derivation it is now clear that the integral 

contribution over Si corresponds to the down-going 
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field at x' . Moreover, for a horizontal free surface 
we know that 0a=O (equation (5)). Since the total 
field consists of a sum of up-going and down-going 
waves, the up-going particle velocity field just below 
5 the surface is therefore given by 

&„(*') = </„(*') + / Vi(x)^ in (x } x')dS(x) . (25) 

Hereinafter the tilde denotes a wavefield 
10 with up-going waves only, i.e. having removed the 
effects of the free surface from the wavefield. 

According to a preferred embodiment of the 
invention, the elastodynamic representation theorem can 
be used to extract the desired up-going wavefield from 
15 surface seismic recordings. For a discussion in the 
case of seismic data from ocean bottom cables, see 
Published UK Patent Application GB 2 333 364 A, 
incorporated herein by reference. 

In the following discussion, 0) is the angular 
20 . frequency, k a =(D/c a and kp = co/cp are the P- and S- 

wavenumbers,. c a and cp are the (isotropic) P- and S- 
velocities, and p is the density. 

The elastic displacement Green's tensor in an 
isotropic homogeneous medium is: 

25 

G,-*(x, x 7 ,^ = -^(k^g^x,**, u)+9^ r (^(x 3 x y , a/)-&(x, x*»)) , (26) 



30 



Where g a and gp are the P- and S-Green's 
functions, and S' m is Kronecker ! s delta. For 
homogeneous media in 3-D these are given by: 
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1 e ik\x-x'\ 

= (27) 



The corresponding isotropic stress Green's 
5 tensor is given in terms of the displacement Green's 
tensor as : 

X ijk - XSi&Gto + + diG jn ) . (28) 

10 A flat horizontal free surface is assumed 

here. Expressions are therefore needed for Green's 
functions such that x = (^,^) and x' =(£',*' 3) are located 
very close to each other in the ^-direction, whereas 
the horizontal separation can be arbitrary. Hence, X3 - 

15 x'3= e -> 0 + and the Green's tensors are shift invariant 
with respect to the X\ and X2 directions. 

Since the "free-space" Green 1 s tensors Gui and 
Z/jjk in equations (26) and (28) are shift invariant, the 
expressions for the up-going wavef ield given by the 

20 representation theorem (25) can be written as spatial 
convolutions denoted by *: 



Dj(x) = I^W-JJWt^x)), (29) 
v 2 (x) = ^(v 2 (x)-F^x)*v i (x))\ (30) 



25 



are: 



The filters in equations (29) , (30) , and (31) 
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= 25S>b(*a = *l) = -24((l + 2^Vc)&+ 2 ^*)- 
J^ = 2S 3v3 (s 3 = 4) =2^(-2^ B + (^ + 2^) 5 f /J ) ) (33) 

where the subscripts v and £ denote either 
index 1 or 2 corresponding to horizontal coordinates x\ 
5 and *2- 

Using the Green's function in equation (27), 
we can obtain explicit expressions for the filters in 
equations (32) and (33) that can be implemented 
directly. In the ^-domain these are: 

10 

*J = irc( 1 - a * 1 (Mt+*f ) <)) ■■ (34) 

K> = -m (i - 2 V 2 (¥%■ + W)) • < 3S) 



where fc 3 (ot) = ^kl-k^ is the vertical 

wavenumber for P-waves and k^ -^k^-k^k^ is the 

15 vertical wavenumber for S-waves. 

A preferred technique for P/S separation of 
surface seismic data will now be discussed. We can 
calculate the divergence and curl of the up-going waves 
by taking spatial derivatives of equations (29), (30) 

20 and (31). . Some care must be taken here since this 
involves spatial derivatives of the stress Green's 
tensor E//*. In this fashion, the following expressions 
are obtained: 
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(V . v) - 'jjfrdivi + d i02 ) + F^\x) * vfa) , (36) 



(V X v) a = -5 3 t*+l*^ <38) 
(Vxv) 3 = ^(9iv a -a 2 y a ). (3?) 

Again, using the Green's function in equation 
(27) , we can obtain explicit expressions for the 
5 filters in equations (36), (37) and (38). In the fk- 
domain these are: 



* - ! *« 2fc (.) 



"3 



(40) 
(41) 



e T)i = i ^ffi^ , (42) 



(43) 



JP< V xv) a = . _ j5i(V xv) t (44) 

10 The traditional way for P-wave interpretation 

is to simply look at the recorded v 3 component . This is 
exact for vertically propagating P-waves. As we saw 
above, the divergence on the other hand is zero at the 
free surface. For horizontally propagating P-waves, the 

15 situation is reversed; V3 is zero whereas the 

divergence exactly contains the P-waves. From equations 
(36) and (40) we see that the correct expression 
combines V3 with the expression for divergence in the 
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presence of the free surface so that it holds true for 
all incidence angles. 

Equations (40) to (44) are all functions of 
the material properties c a and cp at the receiver 
5 location. Estimation of these parameters using a full 
wave equation approach with dense surface measurements 
consistent with that taken here are described in UK 
Patent- Application entitled "System and Method for 
Estimating Seismic Material Properties" (UK Patent 

10 Application No. 0003410.8) filed on 15 February 2000. 

An example of a preferred implementation of 
wavefield separation filters will now be discussed. The 
filters derived above can in theory be implemented 
directly which yield expressions that are exact for 

15 homogeneous media with a flat surface. These 

expressions would remove all down-going and evanescent 
wave types including ground-roll . 

However, the filters are slowly decaying and 
contain some complicating factors. High-order factors 

20 of ki and k2 correspond to high-order derivatives in 

the spatial domain. The most severe problems arise from 

the factors k 3 (ct) and k 3 (P> . These terms do not correspond 

s 

to any straightforward implementation in the spatial 
domain. When they occur in the denominator they 

25 introduce a pole in k a and kp respectively. 

One straightforward filter approximation is 
to make Taylor expansions around = 0 in the 
wavenumber domain (factors of k^ correspond to spatial 
derivatives) . The lowest order terms in the Taylor 

30 approximations to equations (40), (41), (42), (43) and 
(44) are: 
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f(*v) ~ ^-i^^-^), (45) 
Fir^ « (46) 



Figure 2, shows a plot of the normalised real 
part of the filter for removing free surface effects 
5 from divergence in equation (40) . The solid line shows 
the exact filter (equation (40)). For wavenumbers 
beyond the pole corresponding to horizontally 
propagating P-waves, the filter is imaginary 
(evanescent modes) • Figure. 2 also shows a plot of the 
10 zeroth-order (dashed line) , and the first-order (dash- 
dotted line) Taylor approximations as given by equation 
(45). 

The application of a preferred wavefield 
separation approach to synthetic data will now be 

15 discussed. A reflectivity code was used to test the 

wavefield separation approach according to a preferred 
embodiment of the invention. See e.g., Kennett, B. L., 
1983 , . Seismic wave propagation in stratified media: 
Cambridge University Press, Cambridge. This was chosen 

20 as opposed to for instance finite differences since up- 
and down-going wavefields can be calculated separately. 
Moreover, the quantities right below an interface can 
be obtained accurately. The output from the 
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reflectivity code are particle velocities and 
divergence of particle displacement. 

Figures 3a and 3b show an example of a 
seismic survey having locally dense receiver groups, 
5 according to the invention. Figure 3a is also the 
model used to test the wavefield separation approach 
according to a preferred embodiment of the invention. 
In the model, a point-source was used to generate 3-D 
synthetics. The source consisted of a 50 Hz Ricker 
10 wavelet and was of explosive type (radiating P-waves 
only) . 

As shown in Figure 3a source 202 is located 
100 m below earth surface 204. The wavefield was 
recorded at 25 m inter-group spacing from 0 m offset to 

15 450 m. Geophone groups 206, 208 and 210 are shown as 
the first three groups in a series of groups, each 
spaced apart by 25 meters.. At each recording location, 
each receiver group comprises four 3C geophones were 
spaced evenly at 0 . 5 m in both horizontal directions. 

20 Also shown in Figure 3a is a reflective subterranean 
surface 212 that is located 150 meters below earth 
surface 204. 

Figure 3b is a plan view showing an example 
of geophone layout within a single receiver group 210. 

25 Receiver group 210 comprises four geophones 250, 252, 
254, and 256 located on the ground surface. Receivers 
250 and 254 are spaced apart about 0.5 meters from each 
other, and receivers 252 and 256 are spaced apart about 
0.5 meters from each other. Also shown in Figure 3b is 

30 a dashed line 260 on which the other receiver groups 
are positioned, for example, receiver groups 206 and 
208. 
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Advantages of moving to smaller inter- 
receiver spacings include increased accuracy in the 
estimated derivatives, and greater validity in the 
. assumption that the material properties are the same in 
5 the vicinity of the receivers. However, an advantage 
of wider spacing is less sensitivity to noise. These 
competing concerns in combination with the wavelength 
of elastic or acoustic waves (or more accurately, the 
projection of the wavef ield on the recording surface) , 

10 should be taken into consideration when designing the 
receiver group spacing. According to a presently 
preferred embodiment, the locally dense receivers are 
spaced about 1 meter apart. However spacing can in 
some situations be larger, for e.g. around 2 meters, or 

15 smaller, e.g. 0.5 meters. According to a preferred 
embodiment the inter-receiver spacing is around 0.25 
meters or less. As mentioned, the optimal spacing of 
the receivers depends upon the wavelength of interest. 
There should be at least two. receivers within the 

20 projection of the shortest wavelength of interest on 
the recording surface. According to a preferred 
embodiment the receivers are spaced apart by a distance 
approximately equal to or less than one-fifth of the 
shortest wavelength of interest. 

25 In the model based on the geometry shown in 

Figures 3a and 3b, the measurements from the receiver 
groups were used to obtain horizontal derivatives of 
the wavef ield at each location. In this example, we 
will test our wavefield separation techniques on 

30 divergence. We expect that curl should display similar 
results. Finally, all sections shown in this example 
have been plotted using the same scale for amplitudes 
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so that amplitudes can be compared directly to each 
other both between traces and different sections. 

Figures 4a and 4b show the divergence as 
calculated directly by the reflectivity code. Figure 4a 
5 shows a section containing both up and down going 
waves. Note that as vrte approach zero offset the 
divergence vanishes. This is because the divergence is 
zero for vertically propagating plane P-waves. Figure. 
4b shows the divergence of the up-going waves only. 

10 This is the desired P-wave section that we wish to 

obtain using our approach for wavefield separation and 
will therefore serve as our reference solution. Notice, 
that the amplitude variation of different events in the 
section are quite different compared to those in the 

15 section in Figure 4a. Also notice the absence of some 
events due to mode conversions from S- to P-waves at 
the free surface. Some numerical artifacts in the 
numerical solution are also visible (e.g., the flat 
event before the first arrivals) . 

20 Traditionally, 3C data have been interpreted 

by assuming that waves propagate vertically near the 
receivers (steep gradients in material properties are 
assumed in the near- surface region) . Hence, P-waves 
show up on the vertical v 3 component whereas S-waves 

25 appear on the horizontal vi and V2 components. Figures 
5a-d illustrate sections obtained from V3 measurements 
only. To be able to compare them with divergence we 
have applied a time derivative to the V3 measurements 
and scaled them with the P-velocity. Figure 5a shows v 3 

30 divided by a factor of two, according to a conventional 
technique. This exactly corresponds to the up-going P- 
waves for normal incidence (equation (36)). Figure 5b 
shows the difference between this section and the 
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reference solution. We see that the correspondence to 
P-waves rapidly breaks down away from normal incidence, 
where S-waves and mode conversions significantly 
contaminate the result, 
5 Part of the problem is of course that the v 3 

measurements contain both up- and down going waves. 
According to the invention, these can be removed using 
a preferred wavefield separation technique as described 
herein. As can be seen from equations (31) and (35), 

10 this requires knowledge about P- and S-velocities in 

the near surface as well as convolving spatial filters 
with the measured vi and V2 components and adding them 
to the V3 measurements. Figure 5c shows such a 
section. Figure 5d shows the difference between this 

15 solution and the reference solution. Even though the 
result has improved somewhat compared to just taking 
the raw V3 measurements, the result can be improved 
even further. 

Figures 6a-d show the results of using a 

20 wavefield separation technique according to a preferred 
embodiment of the invention. In Figure 6a we have used • 
the zeroth-order Taylor approximation (first term in 
equation (45)). This combines first-order spatial 
derivatives of vi and V2 with a time derivative of V3. 

25 Figure 6b shows the difference between this solution 
and the reference solution. Although, the numerical 
noise that is present in the numerical solution has 
been amplified somewhat by the spatial derivatives, the 
result is now much closer to the reference solution. 

30 In Figure 6c we have used the first-order 

Taylor approximation (both terms in equation (45)). 
This combines first-order spatial derivatives of vi and 
v 2 with a time derivative and second-order spatial 
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derivatives of v 3 . Figure 6d shows the difference 
between this solution and the reference solution. 
Again, the solution has improved even further without 
increasing the noise level compared to the zeroth-order 
5 approximation. 

By comparing all the difference sections in . 
Figures 5 and 6, it is clear that wavefield separation, 
using the techniques according to the invention yield, 
much better results than the traditional P-wave 

10 interpretation techniques. This is particularly the 
case for the window between 0.2 and 0.5 s and 0 m to 
150 m offset in the sections which contain events with 
very realistic incidence angles at moderate to low 
angles with respect to normal incidence (velocity 

15 gradients in the near surface causes energy to be 
incident fairly close to normal incidence) . In this 
example, we estimate that the approximate filters 
improve the results particularly for incidence angles 
up to 30° from normal incidence (this result is 

20 strongly dependent on material properties). However, it 
should be noted that only the theoretical expression 
for the divergence up/down separation filter in 
equation (40) is exact for all wave types. Therefore, 
longer filter approximations will tend to deal better 

25 with horizontally propagating and evanescent wave 
modes. 

According to the invention a new approach for 
P/S separation of land surface seismic data and 
removing the effects of the free surface is provided. 
30 By converting vertical derivatives to horizontal 
derivatives through the use of the free surface 
condition, the methodology can be used even when 
measurements are taken only at the free surface. 
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Therefore by making locally dense measurements of the 
wavefield, all spatial derivatives needed to compute 
divergence and curl of the wavefield at the free 
surface can be obtained. These in turn correspond to P- 
5 and S-waves in isotropic media. 

The effects of the free surface can be 
removed through an up/down separation step as described 
herein. The filter for P-waves depends on both P- and 
S-velocity at the receivers, whereas the S-wave filters 

10 only depend on the S-velocity. Approximations to the 
filters were derived using Taylor approximations and 
tested on synthetic data. 

Filters have been derived for up/down and P/S 
separation using plane wave expansions in Dankbaar, J. 

15 W. M. , 1985, Separation of P- and S-waves : Geophys. 
Prosp., 33, 970 - 986. These expressions can be 
compared to our expressions for full filters. A 
principal difference between the work by Dankbaar and 
those described herein, is that by deploying dense 

20 configurations of 3C geophones, P/S and up/down 

separation in 3D for each recording station can be done 
separately. The preferred up/down separation step makes 
use of approximations to spatial filters to obtain 
operators that are consistent with the number of 

25 geophones in the recording station. This approach is 
more robust since statics and near surface properties 
should be consistent within each recording station. 

The present invention can also be implemented 
in several steps where for instance the P/S separation 

30 is performed in 3D. The up/down separation can be 

performed in 3D or in 2D using an implicit filter if 
data are acquired along a 2D line for instance. 
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For 3C acquisition of land surface seismic 
data it is common practice simply to interpret the 
vertical component as the P-section and the horizontal 
components as SV- and SH-sections. This "traditional" 
5 P/S interpretation is exact for vertical arrivals. 

However, as energy is incident at non-normal incidence 
angles, this approximation tends to break down, both 
because the different waves appear on all components, 
but also because reflection coefficients differ from 

10 unity and mode-conversions occur at the free-surface. 
By comparing the "traditional" P-sections to the new 
methodology using synthetic data, we find a significant 
improvement in obtaining accurate amplitude and phases 
of arrivals for non-normal incidence angles. By simply 

15 using a zeroth-order Taylor approximation, we obtained 
sufficiently accurate results up to incidence angles of 
up to around 30° away from normal incidence (this 
result depends on material properties in the example),. 
Note that a zerbth-order Taylor approximation only 

20 involves first-order derivatives in time and space 
(along the free surface) . Note that the approximate 
expression for divergence consists of two terms. The 
first term corresponds to the divergence in the 
presence of the free surface scaled by a material 

25 constant. The second term is a time derivative of V3 

scaled by a material constant. Hence, a correction is 
added to the "traditional" P-interpretation through the 
first term which improves the accuracy for incidence 
angles outside normal incidence. 

30 Figure 7 is a schematic illustration of a 

seismic data acquisition and processing system, 
according a preferred embodiment of the invention. 
Seismic sources 150, 152, and 154 are depicted which 
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impart vibrations into the earth at its surface 166. 
The vibrations imparted onto the surface 166 travel 
through the earth; this is schematically depicted in 
Figure 7 as arrows 170. The vibrations reflect off of 
5 certain subterranean surfaces, here depicted as surface 
162 and surface 164, and eventually reach and are 
detected by receiver groups 156, 158, and 160. 
Additionally, vibrations in the form of ground roll, 
shown as arrow 128, travel close the surface and can be 

10 recorded by the receiver groups. Each receiver group 
comprises a number of receivers such as in the 
arrangements depicted in Figure 3b. 

Importantly, according to the preferred 
embodiment, the spacing of between the receivers within 

15 a single receiver group is substantially less than the 
spacing between the receiver groups. Schematically, 
this is shown in Figure 7 by the size 122 of receiver 
group 160 is substantially smaller than the distance 
120 between group 160 and an adjacent group 158. In 

20 the example shown in Figure 3a and 3b, the distance 120 
is 25 meters and the size 122 of the receiver groups is 
0.5 meters . 

Referring again to Figure 7, each of the 
receivers in groups- 156, 158, and 160 convert the 

25 vibrations into electrical signals and transmit these 
signals to a central recording unit 170, usually 
located at the local field site. Preferably, the data 
is not group formed, but data from each geophone 
component are recorded. The central recording unit 

30 typically has data processing capability such that it 
can perform a cross-correlation with the source signal 
thereby producing a signal having the recorded 
vibrations compressed into relatively narrow wavelets 
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or pulses. In addition, central recording units may 
provide other processing which may be desirable for a 
particular application. Once the central processing 
unit 170 performs the correlation and other desired 
5 processing, it typically stores the data in the form of 
time-domain traces on a magnetic tape. The data, in 
the form of magnetic tape is later sent for processing 
and analysis to a seismic data processing center, 
typically located in some other geographical location. 

10 The data transfer from the central recording unit 170 
in Figure 7 is depicted as arrow 176 to a data 
processor 180. Data processor 180 can be used to 
perform processing as described in steps 346 through 
352 as shown in Figure 8. 

15 Figure 8 shows steps involved in the 

filtering technique, according to preferred embodiments 
of the invention. In step 340 the seismic survey is 
designed. According to the present invention, the 
survey is preferably designed such that the receivers, 

20 for example geophones or hyrdophones, are positioned in 
locally dense arrangements as shown and described above 
with respect to Figure 3a, 3b, and 7. 

In step 342 a choice is made as to which type 
of wavefield separation is to be carried out. As 

25 described herein, analytical expression can be used for 
up/down separation, p/s separation or both. In the 
case of up/down separation, as described above, one 
would preferably use equations (29), (30), (31), (34) 
and (35) . In the case of p/s separation, as described 

30 above, one would preferably use equations (17) through 
(20) . In situations where both up/down and p/s 
separation is desired, one would preferably use 
equations (36) through (44) . 
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In step 346, an appropriate filter is 
designed based on the desired type of wavefield 
separation (and, accordingly/ the equations selected), 
and based on group recording geometry as designed in 
5 step 340. In creating an appropriate filter one should 
preferably choose a suitable approximation to fit the 
recording geometry. For example, a Taylor series 
approximation is described above, in equations (45)- 
(49). 

10 In step 344, the seismic data is measured 

transmitted, recorded and stored for each individual 
receiver in the survey. 

In step 350, a calculation is made of all 
partial derivatives in time and space needed for the 

15 filter approximations in step 346. 

In step 348, near surface material properties 
are estimated. This can be performed according to a 
number of methods. The presently preferred method is 
using techniques described in UK Patent Application 

20 entitled "System and Method for Estimating Seismic 
Material Properties" (UK Patent Application No. 
0003410.8) filed on 15 February 2000. However, one 
could also interpret near surface properties from the 
surface geology, or alternatively, a shallow refraction 

25 survey, or other suitable means could be used to 
estimate the near surface material properties. 

In step 352 the seismic data is filtered as 
determined in step 346. This is done by combining the 
estimated near surface material properties estimated in 

30 step 348, the seismic data recorded in step 344, and 
the partial derivatives of the wavefield calculated in 
step 350. The result is data 356 which is separated 
as desired (i.e. up/down, p/s or both). 
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While preferred embodiments of the invention, 
have, been described, the descriptions are merely 
illustrative and are not intended to limit the present 
invention. For example, while the preferred 
5 embodiments of the invention have been described 

"primarily for use on the land surface, the invention is 
also applicable to receivers placed on and below the 
ocean floor. In the case of ocean bottom receivers, it 
is preferable to use stress conditions relevant for 

10 fluid- solid boundaries rather than the free surface 
condition. Additionally, . the present invention is 
applicable to seismic measurements made in a borehole, 
known as borehole seismics. Although the examples 
described assume an essentially isotropic medium in the 

15 near surface region, the invention is also applicable 
to anisotropic media. In the case of anisotropic 
media, one may wish to increase the number geophones 
per group. 
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CLAIMS 

What is claimed is: 

5 1. A method of creating a filter for use with 

locally dense seismic data comprising the steps of: 
obtaining survey geometry characteristics 
from a locally dense seismic survey designed to 
record characteristics of an elastic or acoustic 
10 wavefield, the survey comprising a plurality of 

groups of receivers, and each group comprising at 
least three receivers densely spaced from each 
other; 

designing a filter which uses spatial 
15 derivatives of the wavefield of order between or 

including one and the maximum order of spatial 
derivatives of the wavefield that can be estimated 
within a group, such that the filter when 
combining data from within a single group, 
20 separates components of some or all of the 

wavefield arriving at the single group. 

2 . A method according to claim 1 wherein 
the filter is designed to separate up/down going 
25 components of some or all of the wavefield. 

3 . A method according to claim 1 wherein 
the filter is designed to separate p/s components of 
some or all of the wavefield. 



30 



4. A method according to any of claims 1-3 
wherein the filter is designed to separate p/s and 
up/down components of some or all of the wavefield. 
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5. A method according to any of the 
preceding claims wherein each of the densely spaced 
receivers in the group are spaced apart such that 

5 statics in the portion of the wavefield of interest are 
" substantially constant. 

6. A method according to any of the 
preceding claims wherein each of the densely spaced 

10 receivers in the group are spaced apart by about 2 
meters or less. 

7 . A method according to any of the 
preceding claims wherein each of the densely spaced 

15 receivers in the group are spaced apart by a distance 
of about one fifth the shortest wavelength of interest 
or less. 

8 . A method according to any of the 

20 preceding claims wherein the locally dense seismic 

survey comprises generating elastic or acoustic waves, 
and the receivers in a group span less than about the 
smallest wavelength of said elastic or acoustic waves. 

25 9. A method according to claim 1 further 

comprising the step of calculating partial derivatives 
of the wavefield. 



30 



10. A method according to claim 9 wherein 
the step of calculating partial derivatives includes 
using a taylor series expansion as an approximation. 
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11 . A method according to claim 9 further 
comprising the step of estimating near-surface material 
properties . 

12 . . A method according to claim 11 further 
comprising the step of f iltering- the seismic data by 
combining the estimated near surface material 
properties, the seismic data, and the calculated 
partial derivatives. 

13. A method according to claim 1 wherein a 
portion of the wave field that is separated is the 
pressure wavefield. 



15 14. A method according to claim 13 wherein 

the filter separates surface waves and/or air wave 
induced ground motion from the seismic data. 

15. A method according to claim 1 wherein 
20 the locally dense seismic survey is performed on land. 

16. A method according to claim 15 wherein 
the step of creating the filter comprises using the 
free surface condition to convert vertical derivatives 

25 of the wavefield to horizontal derivatives of the 
wavefield 

17 . A method according to claim 1 wherein 
the seismic survey is performed primarily for 

30 hydrocarbon reservoir exploration, evaluation or 
characterisation. 
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18 . A method according to any of the 
preceding claims wherein the near surface velocity is 
essentially isotropic. 

5 19 - A method according to any of the 

preceding claims wherein the near surface velocity is 
anisotropic. . 
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